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Abstract 

A theoretical study is presented of the effect of a radially converging melt flow, which 
is directed away from the solidification front, on the radial solute segregation in 
simple solidification models. We show that the classical Burton-Prim-Slichter (BPS) 
solution describing the effect of a diverging flow on the solute incorporation into 
the solidifying material breaks down for the flows converging along the solidification 
front. The breakdown is caused by a divergence of the integral defining the effective 
boundary layer thickness which is the basic concept of the BPS theory. Although 
such a divergence can formally be avoided by restricting the axial extension of the 
melt to a layer of finite height, radially uniform solute distributions are possible 
only for weak melt flows with an axial velocity away from the solidification front 
comparable to the growth rate. There is a critical melt velocity for each growth 
rate at which the solution passes through a singularity and becomes physically 
inconsistent for stronger melt flows. To resolve these inconsistencies we consider 
a solidification front presented by a disk of finite radius Rq subject to a strong 
converging melt flow and obtain an analytic solution showing that the radial solute 
concentration depends on the radius r as ~ In 1 / 3 (i?o/r) an d ~ ln(i?o/ r ) close to 
the rim and at large distances from it. The logarithmic increase of concentration is 
limited in the vicinity of the symmetry axis by the diffusion becoming effective at 
a distance comparable to the characteristic thickness of the solute boundary layer. 
The converging flow causes a solute pile-up forming a logarithmic concentration 
peak at the symmetry axis which might be an undesirable feature for crystal growth 
processes. 
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1 Introduction 



Solidification and crystallisation processes are present in various natural phe- 
nomena as well as in a large number of material production technologies such 
as, for example, semiconductor crystal growth from the melt, alloy metallurgy, 
etc. Usually the melt used for the production of solid material is not a pure sub- 
stance but rather a solution containing some dissolved dopants or impurities. 
Often the solid material grown from the solution has a non-uniform distri- 
bution of the dissolved substance although the original solution was uniform. 
This non-uniformity is caused by the difference of equilibrium concentrations 
of solute in the liquid and solid phases. Thus, if the equilibrium concentration 
of solute in a crystal is lower than in the melt, only a fraction of solute is 
incorporated from the melt into the growing crystal while the remaining part 
is repelled by the solidification front as it advances into the liquid phase [1]. 
This effect causes axial segregation of the solute, usually concentrated in a 
thin, diffusion-controlled boundary layer adjacent to the solidification front. 
Axial segregation can strongly be influenced by the melt convection. Accord- 
ing to the original work by Burton, Prim and Slichter (BPS) [2], a sufficiently 
strong convection towards the crystallisation front reduces the thickness of 
the segregation boundary layer and so the solute concentration getting into 
the crystal. Such a concept of solute boundary layer has been widely accepted 
to interpret the effect of melt flow on the solute distribution in various crys- 
tal growth configurations [3,4,5]. The BPS approach, originally devised for 
a rotating-disk flow modelling an idealised Czochralski growth configuration, 
supposes the melt to be driven towards the solidification front by a radially 
diverging flow. However, in many cases, as for instance in a flow rotating over 
a disk at rest [6], like in a flow driven by a rotating [7] or a travelling [8] mag- 
netic field, as well as in the natural convection above a concave solidification 
front in the vertical Bridgman growth process [9], the melt is driven away 
from the solidification front in its central part by a radially converging flow. 
Though several extensions of the BPS solution exist (e.g. [10,11,12,13]), the 
possibility of a reversed flow direction away from the crystallisation front has 
not yet been considered in that context. 

In this work, we show that the BPS approach becomes invalid for converging 
flows because the effective boundary layer thickness, which is the basic concept 
of the BPS theory, is defined by an integral diverging for a flow away from the 
solidification front. The divergence can formally be avoided by restricting the 
space occupied by the melt above the solidification front to a layer of finite 
depth, but for higher melt velocities this solution becomes physically incon- 
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sistent, too. Next we consider a solidification front as a disk of finite radius 
immersed in the melt with a strong converging flow and show that a converging 
flow results in a logarithmic solute segregation along the solidification front 
with a peak at the symmetry axis. An analytical solution is obtained by an 
original technique using a Laplace transform. The advantage of this solution is 
its simple analytical form as well as the high accuracy which has been verified 
by comparing with numerical solutions. 

The simulation of dopant transport is an important aspect of crystal growth 
modelling [14,15], and various numerical approaches are used for it. However, 
a numerical approach is always limited in the sense that it provides only 
particular solutions while the basic relations may remain hidden. Besides, 
the numerical solution often requires considerable computer resources when 
a high spatial resolution is necessary which is particularly the case for thin 
solute boundary layers. It has been shown, e.g., by Vartak and Derby [16] that 
an insufficient resolution of the solute boundary layer may lead to numerically 
converged but nevertheless inaccurate results. 

The paper is organised as follows. In Section 2 we discuss the BPS-type ap- 
proach and show its inapplicability to converging flows. The simple model 
problem of radial segregation along a disk of finite radius in a strong converg- 
ing flow is described in Section 3, and an analytical solution for the concen- 
tration distribution on the disk surface is obtained in Section 4. Summary and 
conclusions are presented in Section 5. 



2 Breakdown of BPS-type solutions 

Consider a simple solidification model consisting of a flat radially-unbounded 
solidification front advancing at velocity Vq into a half-space occupied by the 
melt which is a dilute solution characterised by the solute concentration C. 
The latter is assumed to be uniform and equal to sufficiently far away 
from the solidification front. Solute is transported in the melt by both dif- 
fusion with a coefficient D and the melt convection with a velocity field v. 
At the solidification front, supposed to be at the thermodynamic equilibrium, 
the ratio of solute concentrations in the solid and liquid phases is given by 
the equilibrium partition coefficient k. In the absence of convection, the re- 
pelled solute concentrates in a boundary layer with the characteristic thickness 
5 = D/vq. We consider in the following the usual case of a much larger mo- 
mentum boundary layer compared to the solute boundary layer, i.e. a high 
Schmidt number Sc = vjD ^> 1 where v is the kinematic viscosity of the 
melt. The basic assumption of the BPS approach is that the lateral segrega- 
tion is negligible and thus the solute transport is affected only by the normal 
velocity component. The latter is approximated in the solute boundary layer 
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by a power series expansion in the distance z from the solidification front as 
v(z) ~ ^v"(0)z 2 . Then the equation governing the concentration distribution 
in the solute boundary layer may be written in dimensionless form as 

_ 9, dC d?C , . 



where Pe = - ^ <5 ° is the local Peclet number based on the characteristic 
boundary layer thickness So which is used as length scale here while the con- 
centration is scaled by C^. The boundary conditions for the uniformly mixed 
melt core and the solid-liquid interface take the form C\ z ^ (yo — > 1 and 



= 0. (2) 

z=0 



The solution of this problem is 

Pe 



C(z) = 1 + A J exp (-t - -^t 3 ) dt, (3) 



where the constant A = 1 _( 1 J L fc )A(p e ) * s obtained from (2) in terms of A(Pe) = 

J Q °° exp (— t — ^ft 3 ^J dt which according to the relation C"(0) = C ^|p^ < -°- > rep- 
resents an effective dimensionless thickness of the solute boundary layer. Even- 
tually, the concentration at the solidification front is obtained as C(0) = 
[1 — (1 — k)A(Pe)] 1 . This is the central result of the BPS approach stat- 
ing that only the effective thickness of the solute boundary layer defined by 
the local velocity profile is necessary to find the solute concentration at the 
solidification front for a given uniform concentration in the bulk of the melt. 
However, it is important to note that this solution is limited to Pe > and 
it becomes invalid for Pe < when the flow is directed away from the so- 
lidification front because both integrals in Eq. (3) and A(Pe) diverge in this 
case. The goal of this study is to find out what happens to the solute distribu- 
tion when the flow is directed away from the solidification front and the BPS 
solution breaks down. 

The divergence in the BPS model for Pe < is obviously related to the un- 
bounded interval of integration which can be avoided by taking into account 
the finite axial size of the system. The simplest such model, shown in Fig. 1, 
is provided by a flat, radially-unbounded layer between two disks separated 
by a distance 2H. The upper and lower disks represent the melting and so- 
lidification fronts, respectively, and the molten zone proceeds upwards with 
velocity Vq. There is a forced convection in the melt with the axial velocity v(z) 
which is assumed to satisfy impermeability and no-slip boundary conditions. 
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Fig. 1. Sketch of a radially unbounded flat layer with solidification and melting 
fronts at bottom and top, respectively 

There is also a radial velocity component following from the incompressibility 
constraint which, however, is not relevant as long as a radially uniform con- 
centration distribution is considered. Here we choose H as a length scale so 
that the boundaries are at z — ±1. At the upper boundary, there is a con- 
stant solute flux due to the melting of the feed rod with the given uniform 
concentration Cq with velocity vq 



Note that this boundary condition following from the mass conservation does 
not formally satisfy the local thermodynamic equilibrium relating the solute 
concentrations in the solid and liquid phases. In order to ensure equilibrium 
concentrations at the melting front it would be necessary to take into account 
also the diffusion in the solid phase which, however, is neglected here. Such 
an approximation is justified by the smallness of the corresponding diffusion 
coefficient. 

At the lower boundary, coinciding with the moving solidification front, the 
boundary condition is 



Pe (C - Co) + 



dC 

dz 



z=l 



0. 



(l-k)Pe C + 



dC 

dz 



z=-l 



(4) 
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Fig. 2. Modified effective boundary layer thickness PeoA(Peo, Pei) — 1 at the 
solidification front for a horizontal liquid layer of finite height with the flow away 
from the solidification front versus the Peclet number Pe\ of melt stirring at various 
Peclet numbers Peo based on the solidification rate. 



where Pe = v H/ D is the Peclet number based on the solidification velocity. 
The radially uniform concentration distribution depending only on the axial 
coordinate z is governed by 

/ ~ / \ \ dC d C , . 

(-ft. + - = jj, (5) 



where Pe\ is the Peclet number of convection. The solution of the above 
equation is 



Z L 

C(z) = A + B J exp -Pe (t + 1) + Pe l J v{r) dr 



dt. 



(6) 



The boundary condition (4) yields B = —A(l — k)Pe while the remaining 
unknown constant A is determined from the condition at the upper boundary. 
However, for our purposes it is sufficient to express A in terms of the concen- 
tration at the solidification front: A = C{— 1). Then Eq. (6) allows us to relate 
the concentrations at the melting and solidification fronts 

C(-l) = C(l) [1 - (1 - A;)Pe A(P eo , Pe^ 1 , (7) 
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is the effective solute boundary layer thickness defined by the relation 

C A(Pe^Pel) f°H° wm g from Eqs. (4) and (7). This effective boundary layer thick- 
ness at the solidification front is plotted in Fig. 2 for a model velocity distri- 
bution v(z) = (1 — z 2 ) 2 . The effective boundary layer thickness increases with 
the convection but the increase is relatively weak until Pe± becomes com- 
parable to Pe . At this point, the effective thickness starts to grow nearly 
exponentially. 

Although the effective boundary layer thickness is now bounded for any finite 
value of Pei regardless of its sign defining the flow direction, the obtained solu- 
tion is not really free from singularities. At first, note that the concentration at 
the solid-liquid interface becomes singular when the solute boundary layer be- 
comes as thick as -Pe A(Pe , Pe\) = (1 — k)^ 1 resulting in a zero denominator 
in Eq. (7). Second, for larger Pe 1 the denominator in Eq. (7) becomes negative 
implying a negative concentration at the solidification front that presents an 
obvious physical inconsistency. Thus, the obtained solution is applicable only 
for sufficiently weak converging flows and breaks down as the velocity of the 
melt flow away from the solidification front becomes comparable to the growth 
rate at Pe\ ~ Pe . 



3 A disk of finite radius with a strong converging flow 



The assumption underlying both the classical BPS approach and that of the 
previous section is the neglected radial segregation. The simplest physical 
model which could account for radial segregation is presented by a solidifica- 
tion front in the form of a disk of finite radius i? with the melt occupying the 
half-space above it, as shown in Fig. 3. For simplicity, the velocity distribution 
in the melt is assumed to be that of a liquid rotating above a disk at rest. In 
this case, contrary to the classical BPS problem of a rotating disk, the flow is 
radially converging rather than diverging. Thus, within the solute boundary 
layer, assumed as usual to be thin relative to the momentum boundary layer, 
the radial and axial velocity components can be approximated as 
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Fig. 3. Sketch of the solidification front presented by a disk of radius Rq in a 
converging flow. 

Here we choose the thickness of the solute boundary layer based on the axial 
melt velocity as length scale 

do = (2D/v'M) 1/3 , (9) 



and assume the stirring of the melt to be so strong that the advancement of the 
solidification front with the growth velocity Vq is small compared the character- 
istic melt flow in the solute boundary layer. The last assumption implies that 
the local Peclet number based on the growth rate is small: Pe = vodo/D <C 1. 
Then the problem is defined by a single dimensionless parameter, the dimen- 
sionless radius R = Ro/d = R${2D /w"(0)) ^ 3 > which may be regarded as 
Peclet number based on the external length scale Rq and the internal velocity 
scale f = v"(0)d^,/2. The governing dimensionless equation is 





< 10) 



where the radial diffusion term will be neglected as usual for the boundary 
layer solution to be obtained in the following. Sufficiently far away from the 
solidification front a well-mixed melt is assumed with a uniform dimensionless 
concentration Cq = 1. The boundary condition at the solidification front 



Pe (l-k)C + — 



= 0, 

z=0 



8 



for Pe <C 1 suggests to search for the concentration as 

C « C + Pe (l - k)d, (11) 

where C\ is the deviation of the concentration with a characteristic magnitude 
Pe (l — k) <C 1 from its uniform core value C — 1. Then the boundary 
condition for Ci takes the form ^ = — 1, while C is substituted by C\ 

1 dz 2 =0 ' J 

in Eq. (10) which, compared to the original BPS Eq. (1), has an extra term 
related to the radial advection whereas both the term of axial advection due 
to the solidification speed and the radial diffusion term have been neglected. 
Note that on one hand the radial advection term is indeed important because 
without it we recover the BPS case which was shown above to have no bounded 
solution. On the other hand, for the radial advection term to be significant 
the solute distribution has to be radially nonuniform. However, searching for a 
self-similar solution in the form Ci(r, z) = r a F{zr 13 ) leads only to the radially 
uniform solution with a — f3 — 0. This implies that a possible solution has 
to incorporate the radial length scale R. Additional difficulties with finding 
similarity solutions are caused by the explicit appearance of r in Eq. (10). 
Both these facts suggest the substitution r = — ln(r) that transforms Eq. (10) 
into 

/ dC dC\ d 2 C 

z { z te + w) = a* (12) 

with the radial diffusion term neglected as mentioned above. Since the trans- 
formed equation does not explicitly contain r, C (r, z) being a solution implies 
that C{t — r , z) is also a solution. Consequently, r can be replaced by r — r , 
where r = —\n(R) and thus r = \n(R/r). Note that r = corresponds to 
the rim of the disk while r — > oo to the symmetry axis. 



4 Solution by Laplace transform 



Equation (12) can efficiently be solved by a Laplace transform providing 
asymptotic solutions of the solute distribution along the solidification front 
for both small and large r. The Laplace transform defined as C(s,z) = 
Jo 00 C 1 (r,z)e- ST dT transforms Eq. (12) into 

/ dC ^\ d 2 C 



where s is a complex transformation parameter while the boundary condition 
at the solidification front takes the form: ^f- — —-. A bounded solution 



ft-. 



=o 



s 
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of this problem is C(s,z) = ell (§,§,—) , where U(a,b,x) is the confluent 
hypergeometric function [17]. The constant c is determined from the boundary 
condition at the solidification front as c = 3 2 g 3 r^/llj • ^ ^he solidification front 
we obtain 

3^rw3) / £ .i 

( '> s r(2/3) V3'3 



where 



F{ ^=wh- (13) 



The concentration distribution along the solidification front is then given by 
the inverse Laplace transform 



Ci(r,0) = ^ f e st C(s,0)ds 



b—ioo 



The solution for small r follows from the asymptotic expansion of F(p; a) at 
\p\ 3> 1 that can be presented as 



3' 37 V3A3 



where fj(a) are the asymptotic expansion coefficients of F(p; a) = p~ a J2j=o ~r 
which can be found efficiently by the following approach. We start with the ba- 
sic relation F(p; a) = (l + a/p)F(p+l; a) resulting from (13). The asymptotic 
expansion of both sides of this relation can be presented as 

3=0 r j=o P 



where gj(p; a) = (1 + ap x ) (1 + p 1 ) a J = J2i=o 9j ''^ with the expansion co- 
efficients 



mi ) = < 



1, 1 = 

( -^ l (a + j)^ 1 (j + (l-l)(l-a)),l>0 



I I! 



defined by use of Pochhammer's symbol (p) n = T( ^^ ■ Substituting the above 
expansion back into Eq. (14) and comparing the terms with equal powers 
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28 


120 


390 


960 


383040 
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Table 1 "~ ~~ ' 
First 9 coefficients of the series expansion (16) calculated analytically. 

of p we obtain fj(a) = Y^j=o fi( a )gi,j-i( a ), that due to gi$ = 1 simplifies to 
Y4Z0 fi( a )9i,j-i( a ) — 0. Upon replacing j by j + 1 and taking into account 
9j,i( a ) — —j, the latter relation results in 

1 j ~ 1 

fM) = -■ 2 /i( a )^j+i-i( a ). ( 15 ) 



defining fj(a) recursively for j > 0. In order to apply this recursion we need 
/o(a) which can be shown to be constant and therefore fo(a) = 1 because 
/o(0) = 1. Eventually, we obtain 

C 1 (r,0) = f ^g^ +1/3 . (16) 



where dj = njsyl ■ That means, the radial solute segregation along the 
solidification front at the rim is characterised by the leading term Ci(r, 0) ~ 
r(2/3) ln 1 ^ 3 (-R/r). The first 9 coefficients of the series expansion (16) calculated 
analytically by Mathematica [18] are shown in Table 1. The convergence of 

the obtained power-series solution is limited to r < lim^oo ~ 2.09 

[19]. 

The Laplace transform yields also the asymptotic solution for r 3> 1 deter- 
mined by the singularity of the image at s = where 



that straightforwardly leads to 

d(r,0) wco(ln(i2/r) + Cl ), (17) 

where c = « 1.0651, ci = § (> (1) - V (§)) = In + ^= « 0.8516, 

and VK 2 -) is the Psi (Digamma) function [17]. This solution plotted versus 
t = \n(R/r) in Fig. 4 is seen to match both the numerical and the exact power 
series solution (16) surprisingly well already at r > 1. The numerical solution 
of Eq. (12) is obtained by a Chebyshev collocation method with an algebraic 
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Fig. 4. Solute distribution along the solidification front from the rim versus 
r = \n(R/r) resulting from different approximations in comparison to the numerical 
and exact solutions of Eq. (12). 

mapping to a semi-infinite domain for z and a Crank-Nicolson scheme for r 
[20]. 

Note that the solution (17) describing the solute concentration increasing along 
the solidification front as ~ \n(R/r) is not applicable at the symmetry axis r = 
where it becomes singular. This apparent singularity is due to the neglected 
radial diffusion term in Eq. (10) which, obviously, becomes significant in the 
vicinity of the symmetry axis at distances comparable to the characteristic 
thickness of the solute boundary layer (9) that corresponds to a dimensionless 
radius of r ~ 1. The radial diffusion becoming effective at r < 1 is expected to 
limit the concentration peak at ~ ln(R). The asymptotic solution for the solute 
boundary layer forming around the symmetry axis, which will be published 
elsewhere because of its length and complexity, yields for R ^> 1 the peak 
value of the concentration perturbation at the symmetry axis 



Ci(0,0) « c (\n(R) +ci) - c r 



(18) 



where c r ~ 0.3772. The concentration distribution along the solidification front 
in the vicinity of the symmetry axis is shown in Fig. 5. As seen, the solution 
approaches the finite value (18) at the symmetry axis while the asymptotic 
solution (17) represents a good approximation for r > 2. This solution is ob- 
tained numerically by a Chebyshev collocation method [20] applied to Eq. 
(10) with the asymptotic boundary conditions r^- 3 - = z^ 1 = — c 
supplied by the outer asymptotic solution. This defines the solution in the 
corner region at the symmetry axis up to an arbitrary constant which is de- 
termined by matching with the outer analytic asymptotic solution and yields 
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Fig. 5. Concentration perturbation relative to its peak value (18) along the solidi- 
fication front in the vicinity of the symmetry axis together with the corresponding 
asymptotic solution (17). 

the constant c r appearing in Eq. (18). Note that in the described asymptotic 
approximation the difference Ci(r, 0) — Ci(0, 0) shown in Fig. 5 is a function 
of r only while the dependence on R is contained entirely in Ci(0,0) denned 
by Eq. (18). 



5 Summary and conclusions 

We have analysed the effect of a converging melt flow, which is directed away 
from the solidification front, on the solute distribution in several simple solid- 
ification models. First, it was shown that the classical Burton-Prim-Slichter 
solution based on the local boundary layer approach is not applicable for such 
flows because of the divergence of the integral denning the effective thickness 
of the solute boundary layer. Second, in order to avoid this divergence we con- 
sidered the model of a flat, radially- unbounded layer of melt confined between 
two disks representing melting and solidification fronts. This resulted in a ra- 
dially uniform solute distribution which, however, breaks down as the velocity 
of the melt flow away from the solidification front becomes comparable to 
the growth rate. This suggested that a sufficiently strong radially converging 
melt flow is incompatible with a radially uniform concentration distribution 
and, consequently, radial solute segregation is unavoidable in such flows. Thus, 
we next analysed the radial solute segregation caused by a strong converging 
melt flow over a solidification front modeled by a disk of finite radius Rq. We 
obtained an analytic solution showing that the radial solute concentration at 
the solidification front depends on the cylindrical radius r as ~ In 1 / 3 (Rq/v) 
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and ~ In (Ro/r) close to the rim of the disk and at large distances away from 
it, respectively. It is important to note that these scalings do not imply any 
singularity at the axis r = 0. Instead, the concentration perturbation takes 
the value (18) at the mid-point of the finite radius disk. 

It has to be stressed that the radial segregation according to our analysis is by 
a factor ln(i?o/^o) larger than that suggested by a simple order-of-magnitude 
or dimensional analysis (e. g. Eq. (11)). Thus, for converging flows the concen- 
tration at the solidification front is determined not only by the local velocity 
distribution but also by the ratio of internal and external length scales which 
appear as a logarithmic correction factor to the result of a corresponding scal- 
ing analysis. The main conclusion is that flows converging along the solidifica- 
tion front, conversely to the diverging ones, cause a radial solute segregation 
with a logarithmic concentration peak at the symmetry axis which might be 
an undesirable feature for crystal growth applications. 
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